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SECTION  1 


INTRODUCTION 

In  a  previous  work  [3]  we  developed  a  non-local  thermodynamic 
theory  for  the  thermomechanical  behavior  of  a  granular  continuum  in  terms 
of  the  concepts  of  configurational  temperature  and  entropy.  The  theory  was 
applied  to  the  motion  of  granular  media  under  surface  tractions  and  constant 
configurational  temperature. 

The  non-local  theory  differs  from  its  local  counterpart  in  that  the  set  of 
internal  variables  of  the  local  theory  is  replaced  by  an  internal  variable  field 
whose  domain  is  the  interior  of  the  material  body.  The  evolution  equation 
for  the  internal  variables  of  the  local  theory  is  replaced  with  a  partial 
differential  equation  for  the  internal  field.  This  equation  must  be  solved  in 
the  presence  of  the  prevailing  boundary  conditions. 

In  applying  the  theory  to  a  granular  medium  certain  experimental 
observations  had  to  be  considered.  During  a  loading  history  the  applied  load 
is  supported  by  a  sub-set  of  grains  that  form  a  stationary  skeleton.  This 
phenomenon  is  known  as  "arching".  The  skeleton  is  elastically  deformable 
to  the  extent  of  the  elastic  deformability  of  the  individual  grains.  The 
stationary  particles  are  called  the  elastic  phase. 

The  remaining  particles  are  mobile  and  hence  they  are  called  the 
mobile  phase.  The  mobile  particles  move  through  the  elastic  skeleton  in  a 
dissipative  fashion  and  the  displacement  field  of  the  mobile  phase  is,  in  fact, 
the  internal  variable  field  of  the  non-local  thermodynamic  theory. 

The  field  equation  for  ' the  internal  variable  field  is  a  "diffusion 
equation"  with  the  exception  that  the  time  coordinate  is  replaced  by  an 
intrinsic  time,  very  much  along  the  lines  of  endochronic  plasticity.  The 
solution  of  this  equation  gives  the  displacement  field  of  the  mobile  phase. 
This,  and  the  elastic  displacement  field  constitute  the  total  displacement  field 
of  the  granular  continuum.  The  stress  field  is  thereby  determined. 
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During  unloading,  however,  the  entire  phase  is  stationary  and  the 
response  is  elastic  in  consequence  of  the  deformability  of  the  individual 
grains.  Thus  it  is  apparent  that  "mobility"  i.e.,  the  ratio  of  the  mobile  particles 
to  the  total  number  thereof,  in  a  material  region,  depends  on  the  loading 
history  and  may,  in  fact,  depend  on  the  history  of  the  displacement  field  itself. 
The  manner  in  which  mobility  enters  more  generally  as  a  constitutive 
variable  is  a  subject  for  further  research. 

Summary  of  Work. 

In  Section  2  we  lay  the  physical  foundations  of  the  theory  and  develop 
the  concepts  of  non-local  thermodynamics.  Appropriate  field  equations  for 
the  displacement  fields  of  the  elastic  and  mobile  phases  are  derived  together 
with  appertaining  boundary  conditions  necessary  for  the  solution  of  the  field 
equations.  Particle  mobility  appears  implicitly  in  the  constitutive 
formulation  but  is  assumed  uniform  in  the  material  domain  in  question.  It 
does,  however,  vary  with  the  loading  history.  In  a  more  general  formulation, 
mobility  itself  may  be  a  field  with  an  appropriate  field  equation.  Such  a 
theory  will  be  the  basis  of  future  work. 

Section  3  is  devoted  to  the  application  of  the  theory  to  a  granular 
domain  in  the  form  of  a  cuboid  in  compression  with  its  base  rigidly  fixed 
under  conditions  of  uniaxial  strain.  The  resulting  diffusion  equation  is 
solved  under  conditions  whereby  the  intrinsic  time  is  proportional  to  the 
displacement  of  the  top  (mobile)  surface  of  the  cuboid. 

The  mobile  displacement  field  is  thus  determined  at  all  points  in  the 
cuboid.  Comparison  with  experiment  is  made  on  the  plane  near  the  halfway 
point  on  the  height  coordinate  of  the  cuboid.  Agreement  with  experiment  is 
favorable. 

The  solution  was  determined  by  two  different  techniques.  In  the  first 
instance  the  diffusion  equation  was  solved  by  a  finite  difference  method 
using  first  three  and  then  five  equidistant  nodes  along  the  height  of  the 
cuboid.  In  the  second  instance,  an  exact  solution  was  found  by  assuming  the 
height  of  the  cuboid  to  be  asymptotically  large.  This  assumption  was  justified 
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on  the  basis  of  the  solution  obtained  by  the  first  method.  It  was  found  that 
the  two  methods  gave  results  which  were  very  close. 

The  analytical  ("exact")  solution  was  then  used  to  determine  the  stress 
as  a  function  of  the  displacement  of  the  top  surface.  The  calculated  result  was 
then  compared  with  experimental  data  obtained  by  Gill  [4,5]  for  three 
sequential  loading  and  unloading  histories.  Favorable  agreement  with 
experiment  was  obtained. 

In  Section  4  we  deal  with  the  determination  of  material  functions  and 
constants.  Experiments  are  defined  whereby  these  may  be  determined  in  one¬ 
dimensional  as  well  as  three-dimensional  cases,  on  the  basis  of -physically 
realistic  assumptions.  The  details  of  this  method  are  discussed  in  the  body  of 
the  work. 

In  Section  5  we  give  a  conclusive  discussion  of  the  results  and  offer 
recommendations  for  future  research  work. 
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SECTION  2 


THEORY  OF  GRANULAR  MOTION 


KINEMATICS 


The  motion  of  particles  in  a  granular  domain  is  exceedingly  complex. 
Computer  visualization  of  grain  motion  [5]  as  well  as  the  experimental 
observations  of  Gill  [6],  and  Mogami  [7],  indicate  that  there  is  a  subset  of 
particles  that  move  little,  if  at  all,  in  the  course  of  deformation  of  a  confined 
granular  domain.  This  set  of  particles  forms  a  stationary  elastic  skeleton, 
relative  to  which  the  remainder  of  the  particles  move  in  the  course  of  the 
macroscopic  deformation.  It  is  this  motion  that  is  responsible  for  the 
dissipation  associated  with  the  deformation  of  granular  domains. 

We  conclude  that,  in  general,  particle  motion  in  a  granular  domain  is 
described  in  terms  of  an  elastic  and  mobile  phase.  The  elastic  phase  (the 
skeleton)  undergoes  affine  deformation  while  the  motion  of  the  mobile 
phase  is  non-affine  and  is  the  cause  of  energy  dissipation  in  the  domain.  A 
set  of  particles  is  said  to  undergo  affine  motion  if  the  "shell"  formed  by  the 
immediate  neighbors  of  a  particle  is  never  traversed  by  the  particle  during  the 
motion.  Thus,  the  deformation  process  of  an  aggregate  neighborhood  may  be 
thought  to  consist  of  the  deformation  of  the  elastic  "skeleton"  through  which 
mobile  particles  flow  in  a  diffusive  fashion.  This  view  is  in  accord  with  the 
observed  phenomenon  of  "arching"  i.e.,  the  formation  of  skeletal  structures 
that  support  the  applied  stress  without  inhibiting  the  flow  of  the  granular 
medium  as  a  whole. 

It  is  observed  in  aggregates  under  compression,  that  the  overall 
deformation  of  the  domain  is  not  a  continuous  nor  a  homogeneous  process, 
but  consists  of  stochastic  motion  of  individual  particles.  At  any  time  only  a 
subset  of  particles  partakes  in  the  aggregate  motion.  Thus,  excluding  the 
(elastic)  deformation  of  individual  grains,  at  any  time,  the  mobile  particles 
move  relative  to  the  bulk  of  the  aggregate  which  is  essentially  at  rest. 
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In  a  statistical  average  sense,  a  mobile  particle  moves  through  a 
frictional  milieu.  It  is  the  resistance  to  this  relative  motion  that  gives  rise  to 
dissipation  in  the  course  of  the  deformation  of  aggregate  media.  Therefore,  it 
appears  that  the  modeling  of  aggregate  behavior  hinges  on  a  realistic 
description  of  the  diffusive  motion  of  mobile  particles  through  an  aggregate 
domain. 

To  quantify  the  above  ideas,  consider  the  motion  of  particles  in  a 
material  neighborhood,  ideally  an  infinitesimal  region.  Let  N  be  the  total 
number  of  particles,  n  the  number  of  mobile  particles  and  u^  the  (mean) 

displacement  of  the  mobile  particles  in  a  material  neighborhood.  Also  let  Uj 
be  the  mean  displacement  of  the  neighborhood  and  let  the  operator  6  denote 
the  variation  of  a  quantity  it  precedes.  Then: 

5u,  =(n/N)5uf'=5uP  (n/N)  (2.1) 

where  we  have  set  u^  =  uP.  The  "plastic"  connotation  in  5uP  may  not  be  the 
optimal  terminology  because  the  mobile  displacement  is  due  to  translational 
motion  as  well  as  the  deformability  of  the  grains.  However  so  long  as  this  is 
remembered  there  should  be  no  cause  for  confusion. 

Digression.  In  reality  the  stationary  phase  is  elastically  deformable  with 
a  displacement  field  uf.  If  this  field  were  taken  into  account,  then  the 

variation  in  the  mean  displacement  of  the  neighborhood  would  be  given  by 
Eq.  (2.1a): 

<5ui  =  {l  -  (n  /  N)}  5u[  +  (n  /  N)  ^uf'  =  5uj  (2.1a) 


The  most  important  aspect  of  Eq.  (2.1a)  is  that  the  mean  displacement  field  is 
not  a  mere  sum  of  its  elastic  and  plastic  counterparts,  as  in  metal  plasticity, 
but  depends  strongly  on  the  mobility  of  the  grains.  This  is  a  new  and 
important  element  that  must  be  considered  in  the  synthesis  of  constitutive 
equations  for  granular  media. 
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In  the  above  formulation,  account  is  made  of  the  fact  that  n  is  not 
constant  during  deformation.  At  the  present  time,  however,  we  have  limited 
experimental  data  regarding  the  evolution  of  n.  In  fact,  its  measurement  is 
an  interesting  project  for  the  future.  We  do  know,  however,  that  during 
compressive  loading  n  is  of  the  order  of  0.8  N.  This  and  the  fact  that  I  1 6uf  1 1 
«  I  1 6uP  I  I  (where  double  bars  denote  the  norm)  make  the  first  term  on  the 
right-hand  side  of  Eq.  (2.1a)  ignorable  relative  to  the  second,  so  that  during 
compressive  loading: 

Sui  =  5u?  (2.1b) 


A  different  argument  applies  during  unloading,  however.  In  this  case  n  =  0, 
i.e.,  the  particles  are  immobile  and  deformation  of  the  aggregate  is  due  to  the 
elastic  expansion  of  the  the  grains.  Thus  in  the  course  of  unloading: 

5ui  =  (5uf  (2.1c) 

The  displacement  of  the  mobile  particles  relative  to  the  neighborhood 
is  then  Uj  where, 

^uj  =  5uP  -  5uj  =  (1  -  n  /  N)  (2.2) 

As  we  shall  see  Uj  is  important  in  the  determination  of  the  resistive  force  in 
course  of  particle  motion.  It  follows  from  Eqs.  (2.1)  and  (2.2)  that: 

5u-  =  (1  -  n  /  N)(5uP  (2.2a) 

THERMODYNAMICS 

At  present,  irreversible  thermodynamics  deals  with  the 
thermodynamic  state  of  a  neighborhood.  Heredity  is  dealt  with  in  terms  of 
locally  homogeneous  neighborhoods  whose  thermodynamic  state  is  specified 
by  the  (locally  uniform)  strain,  and  a  requisite  number  of  internal  variables 
which  are  also  locally  uniform.  This  approach  leads  to  local  constitutive 
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theories  according  to  which  uniform  boundary  tractions  give  rise  to  uniform 
deformation  gradients.  However,  a  different  and  more  realistic  point  of  view 
is  possible  if  one  regards  the  entire  material  domain  as  a  thermodynamic 
system,  rather  than  attempt  to  infer  its  behavior  from  that  of  an  idealized 
(locally  uniform)  material  subdomain. 

Thermodynamics  of  the  Mobile  Phase 

We  begin  with  a  review  of  thermodynamics  of  internal  variables 
under  isothermal  conditions.  Let  Xj  denote  the  "observable"  kinematic 
variables  and  Xj  the  corresponding  dual  forces.  Also  let  qi  denote  the  internal 
variables  and  Qi  the  dual  internal  forces.  If  M/(xi,qi)  is  the  Helmholtz  free 
energy  of  the  system,  then; 


Xj  =  (^Xjlq 

(2.3a) 

Qi  =-<?(/// ^qj  lx 

(2.3b) 

In  the  case  of  linear  thermodynamics,  following  Biot  [7],  Qi  are  related 
to  the  rates  of  change  of  qi  linearly,  through  an  "Onsager"  positive  definite 
and  symmetric  matrix  bij  so  that, 

Q,  =  bjj  (9qj  /  dz  (2.4) 

where  z  is  an  intrinsic  time  scale,  not  necessarily  Newtonian.  The  concept  of 
intrinsic  time  was  first  used  by  Valanis  [9,10]  in  the  development  of 
Endochronic  plasticity.  Equations  (2.3b)  and  (2.4)  combine  to  give  the 
evolution  equations  for  qi  in  the  form: 

dxf/  /  ^qj  +  bjj  (9qj  /  ^z  =  0  (2.5) 


We  note  that  Eqs.  (2.3a)  and  (2.3b)  give  rise  to  the  relation: 


I  Xj  5xj  -  X  Qi  5qi  =  5\j/ 


(2.6) 
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where  5  denotes  an  arbitrarily  small  variation  in  the  quantities  that  it 
precedes. 

In  the  past,  constitutive  theory  was  developed  on  the  basis  of  local 
thermodynamics  of  a  vanishingly  small  neighborhood  in  which  event  Xi 
denotes  the  components  of  the  strain  tensor  of  the  neighborhood  while  Xj 
denotes  the  components  of  the  dual  stress  tensor  in  the  sense  that, 

XidXi=dW  (2.7) 

where  dW  is  the  increment  of  work  per  unit  volume. 

Global  Thermodynamics  of  the  Mobile  Phase 

In  the  case  of  global  thermodynamics  developed  here,  the  entire 
material  domain  constitutes  the  thermodynamic  system.  The  observable 
variables  are  the  surface  displacements  UP  while  the  internal  variables  are  the 
displacements  uP  in  the  interior  of  the  domain.  The  Helmholtz  free  energy  T 
of  the  system  is  the  integral  of  the  free  energy  density  y/  over  the  material 
domain  D.  In  turn  y/  is  a  function  of  the  gradient  of  the  displacement  field 
uP  just  like  in  the  case  of  the  local  theory. 


Thus,  insofar  as  the  mobile  phase  is  concerned, 

V^=V'(uf,j)  (2.8) 


D 


In  the  case  of  the  present  thermodynamics  system  Eq.  (2.6)  becomes: 

J  Tf  Suf  dS  -  Jq,  5uf  dV  =  (54^  (2.10) 

S  D 
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where  TP  are  the  surface  tractions  on  the  mobile  phase. 


We  remark  that  X  Xj  6xi  in  Eq.  (2.6)  is  an  increment  of  work  done  by  the 
external  forces  so  that  if  a  body  force  field  were  present,  then  Eq.  (2.10)  would 
be  modified  to  include  the  effect  of  such  a  field  in  which  case  it  would  read: 


(2.10a) 


In  the  case  of  the  mobile  phase  the  only  possible  such  field  is  the  gravitational 
field  which  we  ignore.  Thus  fP  =  0. 

The  evolution  equation  for  the  internal  variables  now  becomes: 

Qi=bij^P/<9z  (2.11) 


To  obtain  the  resulting  thermodynamic  equations  we  need  to  derive  an 
explicit  expression  for  S'T.  It  follows  from  Eq.  (2.9)  that, 

54^  =  J(^t///<9uP,j)5uP,j  dV  = 

D 

(2.12) 

j[(d¥/  )^uP},j  dV  +  /  ^P,j)5uP  dV 

D  D 

When  the  Green  Gauss  theorem  is  applied  to  the  first  integral  on  the  right- 
hand  side  of  Eq.  (2.12)  and  use  is  made  of  Eq.  (2.1),  the  following  relation 
results: 
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5ufdW  =  0 


(2.13) 


where  use  was  made  of  the  fact  that  on  the  surface  uP  =  UP.  Since  the 
variations  5UP  and  6uP  are  arbitrary,  the  integrands  in  Eq.  (2.13)  must  vanish 
separately,  in  which  case  the  following  thermodynamic  relations  result: 


Tj  =1^1/// ^P,j|nj  onS 

Qi  inD 


(2.14) 

(2.15) 


Thus  using  Eq.  (2.11)  one  obtains  Eq.  (2.16)  which  is  the  evolution 
equation  for  the  internal  variables,  in  the  form  of  a  partial  differential 
equation  in  the  spatial  coordinates  yj  and  the  intrinsic  time  z: 

/  ^P,j|,j  =  bij  duf  /  dz  .  (2.16) 


In  the  event  that  the  material  is  isotropic  (in  a  continuum  sense)  then  b  is  an 

isotropic  second  order  tensor,  and  in  the  case  where  it  is  a  constant,  it  is 
proportional  to  the  unit  matrix.  In  this  study  we  take  the  position  that  b  is 

given  by  Eq.  (2.17). 

'  (2.17) 


where  b  is  a  scalar  invariant.  Equation  (2.16)  becomes, 
|<9V'  /  /j|,j  =  b  ^P  /  dz 


(2.18) 


Equation  (2.18)  is  the  evolution  equation  for  the  internal  variable  field  uP.  It 
is  a  field  (diffusion)  equation  which  must  be  solved  in  the  light  of  the 
prevailing  displacement  or  stress  boundary  conditions  and  as  such  it  is  a 
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GLOBAL  equation  Note  that  in  local  thermodynamics  the  evolution 
equation  for  the  internal  variables  is  independent  of  the  boundary  conditions 
of  the  material  domain. 

The  Resistance  Coefficient  b.  Physically,  resistance  is  generated  by 
differential  particle  motion,  in  our  case  defined  by  uf  in  Eq.  (2.2).  Thus  the 
resistance  force  Qi  is  given  by  Eq.  (2.19). 


Qi  =  b  ^*  /  <9z  (2.19) 

where  b  is  a  material  property  which  depends  on  the  mobility  of  the  grains 
n/N,  their  geometry  and  the  nature  of  the  grain  surface  but  also  on  the 
confining  pressure. 

In  the  non-local  theory  the  lateral  stress,  brought  about  by  lateral 
confinement,  varies  with  the  height  of  the  cuboid.  This  poses  difficulties  in 
solving  the  prevailing  equations  analytically.  Analytical  solutions  are 
necessary,  if  the  experimental  data  are  to  be  utilized  to  obtain  the  material 
constants  and  functions.  To  this  end  we  have  set  b  to  depend  on  the  mean 
pressure  of  the  entire  domain  of  the  cuboid,  and  therefore  on  the  applied 
stress  o.  In  addition  b  will  depend  on  the  porosity  w  and  therefore  on  the 
volumetric  strain.  This,  in  the  uniaxial  case,  is  one  third  of  the  axial  strain. 

Remark.  Thought  was  given  to  the  question  of  how  porosity  depends 
on  the  strain  of  a  granular  neighborhood.  Evidently  it  depends  on  the 
positional  configuration  of  the  grains  as  well  as  their  deformation.  For,  if  the 
external  boundary  of  granular  domain  were  held  fixed  and  the  grains  were 
compressed  elastically,  the  porosity  would  obviously  increase.  In  a  loading¬ 
unloading  sequence  the  porosity  would,  therefore,  depend  on  the  cumulative 
strain  of  the  neighborhood. 

Because  pressure  and  porosity  dependence  are  due  to  two  independent 
physical  processes,  b  will  depend  on  these  in  a  multiplicative  fashion.  Hence: 

b  =  b'  f(a)  g(w)  (2.20) 
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The  constant  b’  is  a  property  of  the  granular  medium.  We  now  reason  that  Qj 
is  proportional  to  the  number  of  friction  surfaces  and,  therefore,  proportional 
to  n.  Hence: 


b'  =  nbo  (2.20a) 

Thus  in  view  of  Eqs.  (2.19),  (2.20)  and  (2a),  the  dissipative  force  Qi  will  be 
given  by  Eq.  (2.21): 

Q.  =  n{l  -  n  /  N)}bo  f{a)  g(w)  9uf  /  dz  (2.21) 

where  f  and  g  are  material  functions  to  be  determined  by  experiment. 

The  Stress  Distribution  in  a  Granular  Domain. 

Attention  must  be  given  to  this  aspect  of  constitutive  behavior  of  a 
granular  medium.  There  are  in  fact  many  subtle  points  which  are  not 
immediately  evident  upon  initial  examination.  The  essence  of  the  problem 
lies  in  the  fact  that  the  stress  is  transferred  from  the  mobile  phase  to  the 
elastic  skeleton  by  friction.  It  is  the  friction  that  inhibits  the  collapse  of  a 
granular  aggregate  when  the  latter  is  subjected  to  surface  tractions. 

The  simplest  possible  model  that  bears  resemblance  to  the  physics  of 
the  motion  is  that  of  a  friction  bjock  (the  mobile  phase)  moving  between  two 
elastic  walls,  which  represent  the  elastic  skeleton.  See  Figure  1.  The  block  is 
of  much  larger  cross-sectional  area  than  that  of  the  walls,  since  the  mobile 
particles  are  much  more  numerous  than  their  stationary  counterparts.  The 
depth  of  the  block  reflects  the  extent  to  which  the  stress  has  been  transmitted 
to  the  interior  of  the  mobile  phase. 

One  can  see  that,  at  the  surface,  most  of  the  stress  is  carried  by  the  block, 
(the  mobile  phase).  However,  in  the  interior,  and  in  a  manner  consistent 
with  experimental  data,  the  skeleton  carries  a  larger  stress  commensurate 
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with  the  depth,  as  stress  is  transferred  from  the  block  to  the  elastic  skeleton 
(the  walls)  by  friction. 

Thus  letting  Tf  denote  the  surface  traction  on  the  elastic  phase,  TP  the 
traction  on  the  mobile  phase,  and  Ti  the  traction  on  the  surface  of  the 
granular  medium,  then  by  simple  proportion; 

T;^  =  {l-(n/N)}Ti  ;  T|P=(n/N)Ti  (2.22a,b) 


Hence; 

Ti=Tf  +  T|P  (2.23) 

Some  further  obser\'ations  are  in  order.  It  follows  from  Eq.  (2.14)  that; 

=  (2.24) 

where  oj-  is  the  stress  in  the  mobile  phase.  Thus  Eq.  (2.15)  reads; 

CTP,j  =  Q,  (2.25) 


Thus  Qi  acts  as  a  negative  body  force  brought  about  by  the  interaction  between 
the  mobile  and  elastic  phase.  Since  action  and  reaction  are  equal  and 
opposite,  an  equal  and  opposite  force  must  be  exerted  on  the  elastic  phase.  It 
follows  that  the  equilibrium  equation  for  the  elastic  phase  must  read; 

C75,j+Qi=0  (2.26) 


Furthermore  as  a  result  of  Eqs.  (2.25)  and  (2.26); 


Cr®  +  (T?  =  <7;; 
1)  1)  4 


(2.27) 
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(2.28) 


CTij,j  =  0 


where  aij  is  the  resultant  stress  in  the  granular  medium. 
THERMODYNAMICS  OF  THE  ELASTIC  PHASE 

We  conclude  this  subsection  with  the  treatment  of  the 
thermodynamics  of  the  elastic  phase.  Grain  mobility  and  deformability  are 
independent  processes  in  the  sense  that  they  pertain  to  disjoint  subsets  of 
grains.  It  is  reasonable,  therefore,  to  expect  the  free  energy  density  \\i  of  the 
®8§^®8^te  domain  to  be  the  sum  of  the  individual  densities  of  the  elastic  and 
mobile  phase.  Thus: 

¥=¥p  +  ¥e  (2.29) 


where 


V^p  =  ¥p 


(2.30) 


while 


¥e  =  ¥e 


)  J 


(2.31) 


The  essential  difference  between  the  elastic  and  the  mobile  phase  is  that  in 
the  elastic  phase,  the  dissipation  is  zero,  while  a  body  force  field  Qi  arises  as  a 
result  of  the  dissipation  in  the  plastic  phase. 

Thus  in  the  case  of  the  elastic  phase  the  variational  equation  (2.10a) 
will  read 
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(2.32) 


since  in  fact  fi  =  Qi.  Following  the  technique  previously  applied  to  the  plastic 
phase  we  find  that 

nj  on  S  (2.33) 


dy/ 


dy/ 


+Qj  =0  in  D 


(2.34) 


of  course  we  recognize  0vi//3uj®  as  ofj  so  that  Eq.  (2.34)  is  in  conformance  with 
Eq.  (2.26)  which  was  obtained  by  static  considerations. 


The  Displacement  Field 


Following  the  previous  notation  uf  is  elastic  displacement  field  pertaining  to 
the  skeleton,  where  uP  is  the  displacement  field  of  the  mobile  phase.  Let  Ui  be 
the  displacement  field  of  the  aggregate.  Evidently,  at  the  surface,  the  three 
fields  must  satisfy  the  condition 


(2.35) 


at  this  point  we  utilize  Eqs.  (2.22a,b)  to  obtain 
(l-n/N)5u^ +^5uP  =  5ui 


(2.36) 


which  is  identical  to  Eq.  (2.1a).  Equation  (2.36)  relates  the  displacement  field 
of  the  aggregate  to  the  displacements  of  the  individual  phases. 
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It  is  observed  that  the  skeleton  is  not  static.  It  changes  configuration 
and  is  not  always  made  up  of  the  same  particles.  Nonetheless  on  physical 
grounds  the  compliance  of  the  skeleton  is  directly  related  to  the  deformability 
of  the  individual  particles.  The  deformation  of  the  skeleton  is  small 
compared  to  the  aggregate  deformation  due  to  particle  motion.  Therefore  we 
expect  that 


(2.37) 


and  since,  furthermore,  the  number  of  mobile  particles  n  is  much  larger  than 
that  of  the  stationary  counterparts  (as  estimated  values  of  n/N  is  0.8),  one 
may  ignore  the  term  (1  -  n/N)uf  relative  to  n/N  u?  in  Eqs.  (2.36).  Hence,  to  a 
good  degree  of  accuracy 


(2.38) 


Equation  (2.18)  in  One  Dimension 


In  one  dimension  (uf  =  uP,  uf  =  uf  =  0)  Eq.  (2.18)  becomes: 


^(jP  _  ^  (9uP 
(?x  Bz 


(2.39) 


where 


£rP= _ _ 

3(auP/<?x| 


Within  the  restriction  of  small  strains 


A 


2 


(2.40) 


(2.41) 
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where  A  is  an  appropriate  elastic  stiffness  coefficient.  (In  general  A  will 
depend  on  the  confining  pressure  as  well  as  on  the  porosity  and,  therefore,  on 
the  plastic  volumetric  strain  (3uP/3x)).  It  follows  from  Eqs.  (2.39)  and  (2.40) 
that: 


(jP  =  A 


1  ^A 

2  (9^^P  / 


(2.42) 


Since  in  the  linear  theory  terms  in  (3uP/3x)^  are  ignorable,  then 

(7P  =  A-^  (2.43) 

To  make  matters  simple,  the  dependence  of  A  on  a  and  w  is  assumed  to  be  of 
the  same  as  that  of  b.  Thus: 

A  =  Aof(a)g(w)  (2.44) 

and  hence: 


ctP  =  Aof(cr)g(w)^  (2.45) 

At  this  point,  we  substitute  Eq.  (2.45)  in  Eq.  (2.39),  retain  linear  terms  and  use 
Eq.  (2.21)  to  find: 

A„-^  =  N(l-n/N)b„^  (2.46) 

dx 

This  is  the  equation  that  we  shall  use  in  the  determination  of  the  internal 
displacement  field  in  simple  compression  generated  experimentally  by  Gill 
[4,6]. 
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Explicit  Three-Dimensional  Equations 


In  the  context  of  small  strains  and  an  initially  isotropic  granular 
medium  we  set 

(.2.47) 

where  the  elastic  constants  X  and  p  will  depend  on  the  porosity  as  well  as  the 
prevailing  hydrostatic  stress.  We  now  substitute  Eq.  (2.47)  in  Eq.  (2.48), 
perform  the  indicated  differentiation  and  retain  only  the  linear  terms  in  the 
displacement  gradient  to  obtain  the  following  equations  for  the  plastic 
displacement  field  uP 


(A+p)uP,ji+puP,jj  =  b 


~dr 


(2.48) 


We  now  express  uP  in  terms  of  its  (independent)  irrotational,  and 
isotropic  components  Ui  and  Uj*,  i.e., 

uP  =  ui  +  u*  (2.49) 

such  that 

Uj  ,j  —  0  ,  Giji<  ^k,j  ~  ^  (2.50a,b) 


where  eijk  is  the  permutation  symbol.  This  decomposition  allows  Eq.  (2.48)  to 
be  expressed  into  two  separate  relations.  Using  Eq.  (2.50a, b)  Eq.  (2.48)  becomes: 


(A+p)ui,ji 


+  b 


du* 

dz 


in  view  of  Eq.  (2.50a).  However,  in  view  of  Eq.  (2.50b) 


(2.51) 
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Thus,  Eq.  (2.51)  becomes 


Thus  the  motion  is  decomposed  into  a  deviator  field  and  a  volumetric  field. 
The  deviator  field  obeys  Eq.  (2.54)  while  the  isotropic  field  obeys  Eq.  (2.57). 

Probabilistic  Considerations 

In  support  of  the  thermodynamic  arguments  discussed  previously  we 
pursue  a  probabilistic  theory  akin  to  the  one  that  governs  diffusion  in  liquids 
with  Brownian  motion  [11].  To  this  end  let  f(xi)  be  a  property  of  the  medium 


where  Xj  are  the  coordinates  of  the  particles  of  the  medium.  Let  a  disturbance 
cause  a  change  in  Xi  which  we  denote  by  Axj,  in  time  At.  We  note  that 
microscopic  changes  are  discreet  and  that  At  is  the  average  minimum  time 
over  which  a  microscopic  change  can  take  place.  Thus, 

Axi=ViAt  (2.58) 


where  At  s  i  is  a  material  property  of  the  process,  a  "characteristic  time,"  so  to 
speak,  and  Vi  is  the  particle  velocity  vector.  The  field  Axj  is  not  deterministic, 
(i.e.,  it  is  random)  in  the  sense  that  the  probability  of  a  velocity  Vj  is  p,  where 

P  =  </>(v,)  •  (2.59) 

We  now  assume  that  motion  in  three  mutually  perpendicular  directions 
constitutes  three  independent  events,  in  which  case, 

</>(Vi)  =  (/>(vi)(^(v2)0(v3)  (2.60) 


and  let 


(p{y)dv  =  1 


CO 

J* V  0(v)dv  =  V 

— oo 


oo 

J*  v^0(v)dv  =  v^ 

— oo 


(2.61a,b,c) 


In  the  absence  of  externally  applied  body  forces,  we  stipulate  that  (l)(v)  is 
a  symmetric  function  in  the  sense  that  forward  and  backward  motion  are 
equally  likely.  Thus, 


v  =  0 


(2.62) 


We  note  that  v^  is  a  property  of  the  distribution  (j)(v). 

Following  a  disturbance,  the  expected  value  f  of  f  is: 
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eo 

Hi 


f=  I  f(xi+ViAt)0(v,)dvi 


(2.63) 


Note  that  the  expected  value  fo  iri  the  undisturbed  configuration  is 


oo 

fo=  JJ  Jf(xi)0(vi)dVi  =f(xi; 


(2.64) 


in  view  of  Eq  (2.61a).  Thus, 


A?  =  f  -  fo  =  J  J  J  {f(xi  +  VjAt)  -  f(xi  )}(/>(vj  )dvi 


(2.65) 


Limiting  the  study  to  disturbances  that  are  of  small  order,  say  0(e): 


1 

f(xi  +  Au j )  -  f(xj )  =  f , j  Au j  +  —  f ,ij  Au j  Auj  +  0(e  ) 


(2.66) 


Thus, 
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oo 


df 


v2  +  rp-V3  0(v2)0(v2)^(v3)dvidv2dv3 
^3  J 


V2  + 


A 

0(vi)0(v2)(^(v3)dvidv2dv3 

J 


V1V2  + 


^1^X3 


-V1V3  + 


Bx2^3 


V2V3 


</’(vi)0(v2)0(v3)dVidV2dV3 


which  in  view  of  Eq.  (2.72),  becomes. 


£2 

2! 


^  d^i 


d^i'^ 


i?X?  <?Xo 


3; 


v‘ 


But 


7 

Af  =  — At  =  — T 

dt  dt 


(2.67) 


(2.68) 


(2.69) 


and  hence: 

■^  =  DV2/  (2.70) 

where  denotes  the  Laplacian  operator  and 
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(2.71) 


i.e.,  the  diffusion  coefficient. 

However,  t  is  a  characteristic  material  time.  Thus  we  may  define  an 
intrinsic  time  z,  such  that 


T 


in  which  case  Eq.  (2.70)  becomes: 


dz 


a  T^f 


(2.72) 


(2.73) 


where 


a  = 


(2.74) 


Specifically  in  the  case  where  f  is  the  isochronic  displacement  fields  uf,  then 


dz 


=  a  V^uj 


(2.75) 


where 


(2.76) 


while  if  ?  is  the  irrotational  displacement  field  potential  <{),  then 


dz 


(2.77) 
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where 


(2.78) 

b 

Thus  probability  theory  and  irreversible  thermodynamics  have  been  brought 
into  complete  correspondence. 
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ELASTIC  SKELETON 


Applied  Stress 


■ 

II 

mn 

1 

1 

1 

■ 

■ 

FRICTION  BLOCK 


STRESS'  FREE 
REGION 


Fig.  1.  Friction  Block  Model 
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ELASTIC  SKELETON 


SECTION  3 


APPLICATIONS 


DISPLACEMENT  FIELD  UNDER 
CONSTRAINED  AXIAL  COMPRESSION 


Case  (1)  Monotonic  Loading. 


The  theory  has  been  applied  to  the  case  of  vertical  motion  of  particles 
in  a  cuboid  under  conditions  of  axial  strain.  Insofar  as  experiment  does  not 
predict  uniform  strain  under  uniform  compressive  stress,  the  theoretical 
consequences  were  put  to  the  test  by  solving,  with  a  finite  difference 
technique  the  diffusion  equation  derived  from  Eq.  (18),  for  the  particle 
displacement  at  the  midway  plane  of  the  cuboid.  Excellent  agreement  was 
obtained  between  theory  and  measurement.  This  outcome  is  encouraging  in 
that  it  supports  the  conceptual  framework  of  the  theory  which  is  essential  to 
the  understanding  of  the  deformation  of  large  bulks  of  granular  domains. 

The  same  technique  was  also  used  to  obtain  theoretically  the  particle 
displacement  at  points  located  at  one  quarter  and  three-quarters  of  the  height 
of  the  cuboid.  Experimental  results  at  these  points  were  not  available  for 
comparison.  An  analytical  solution  was  also  obtained  by  an  "infinite  length" 
approximation.  This  solution  is  in  closed  form  and  also  agrees  with 
experiment. 

Finite  Difference  Scheme. 


We  proceed  to  solve  by  a  finite  difference  technique  the  one¬ 
dimensional  generic  diffusion  equation 


d^u  _  ^ 


(3.1) 


where 
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a  =  (Ao/boXl-n/Nr'N-'  (3.1a) 

Let  the  height  h  of  the  cuboid  be  divided  into  N+1  equal  intervals,  with  a 
resulting  number  of  N  interior  nodes  and  two  boundary  nodes  one  at  each 
boundary. 

The  differential  Eq.  (3.1)  is  then  satisfied  at  all  interior  points,  thus 
giving  rise  to  a  set  of  simultaneous  ordinary  differential  equations  shown 
below: 


(2ui-U2)  + 


1  duj 
Cn  dz 


(-ui  +  2u2  -U3)  + 


1  du2 
Cjsj  dz 


=  0 


(3.2) 


(-Un_,  +  2un-U)  +  ^^  =  0 

L-Inj  dz 


where 


[n/(N  +  l)]^ 


Equations  (3.2)  may  be  written  in  the  matrix  form 


(3.3) 


(3.4) 


where 


2 

-1 


[A]  = 


-1 

2  -1 

-1  2 


-1 


2 


and 


0 

0 

6 


There  exists  a  diagonalizing  matrix  Rjj  where 


1 


rw+n 


i/. 


sin 


ij  ;r 

JT+i 


\ 

J 


y  2  J 


such  that 

r'''ar  =  [a] 


A 


n 


=  4  sin 


2 


n  n 

2(N  +  1) 


(3.5) 


(3.6) 


(3.7) 


(3.8) 

(3.9) 


Define: 

u*  =  Ru  ,  u  =  Ru  (3.10a,b) 

Cn  r\u  + —  =  r'^u  =  U* 

^  J  dz  (3.11) 

Then,  in  view  of  Eqs.  (3.8)  and  (3.9)  the  following  set  of  equations  result: 
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(3.12) 


- 

Ul 

dui 

2-2 

* 

U2 

1 

dz 

. 

+ - 

=z 

*• 

cn 

dUn 

dz  J 

which  are  N  uncoupled  ordinary  differential  equations  of  the  form 


Aj.  Uf  + 


1  du* 
Cjsj  dz 


U^Cz) 


(3.13) 


Knowing  U?  one  may  then  solve  Eq.  (3.13)  and  use  Eq.  (3.10b)  to  determine  Ur 
as  a  function  of  z  all  the  modes  r  =  1  ...  N. 

To  calculate  theoretically  the  mean  vertical  displacement  of  particles 
located  points  h/4,  h/2  and  3h/4  from  the  base  of  the  cuboid,  we  subdivided 
the  height  h  into  four  equal  intervals  with  a  resulting  number  of  three 
interior  nodes.  Thus  in  Eq.  (3.3),  N  =  3. 

Case  of  N  =  3 


Use  of  Eqs.  (3.7)  and  (3.8)  gives  rise  to  the  following  results: 


1 

V2 

1 


V2  1 

0  -^^2 

-V2  1 


(3.14) 


[A]  = 


1- 


2/2' 

s  . 


1  + 


V2 


(3.15) 


Also, 
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{h/4f 

Solution  of  Eg.  (3.13) 


(3.16) 


It  follows  from  Eqs.  (3.6),  3.10b)  and  3.14)  that 

U;  =  brU  (3.17) 


where 


br 


:I1  ll 

2  2j 


(3.18) 


In  view  of  our  discussion  in  Section  2,  under  monotonic  conditions, 


d2  =  d 


f— 1 


(3.19) 


In  keeping  with  the  linear  form  of  the  theory  we  take  z  to  be  its  average  value 
z  over  the  entire  domain  of  the  cuboid.  Thus 


dz  =  dz  = 


dx 


dx  = 


dU 


(3.20) 


and  hence: 

dz  =  -^  (3.21) 

h 

at  all  points  along  the  height  of  the  cuboid. 

In  view  of  Eq.  (3.21),  Eq.  (3.13)  becomes:; 
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(3.22) 


Aj.  Uj  + 


_h^ 

c*  dU 


=  brU 


This  has  been  solved  by  Laplace  transform  in  the  light  of  the  initial  condition 
uf(0)  =  0,  to  give 


(3.23) 


where 


We  remark  that  the  only  material  constant  in  the  above  equation  is  c*. 
Evidently  Ur  may  be  obtained  from  u*  by  the  relation 
Ur  =  Rrs  U^  (3.25) 

following  Eq.  (3.10b)  and  the  fact  that  R  = 

Since  Eq.  (3.1)  is  a  generic  equation  we  obtain  the  solution  to  our 
problem  by  setting 

u^  =  Uj.  (3.26) 

Computations 

It  is  clear  from  Eq.  (3.33)  that  the  only  material  constant  is  c*.  To  obtain 
c*  we  made  use  of  experimental  measurements  made  by  Gill  [4]  at  a  "point" 
close  to  the  mid-height  of  the  cuboid,  and  derived  the  theoretical  expression 
for  N  =  1.  Omitting  the  analysis  which  was  given  previously  [2], 
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(3.27) 


where  U  is  the  measured  displacement  and  c  =  1  /2  Specifically 


c  =  8^^ 


bolN-n; 


(3.28) 


Good  agreement  with  measured  values  of  uP  was  found  when  c  was  set  equal 
to  0.6,  i.e.,  c*  =  1.2. 

When  uP  was  calculated  at  h/2  for  this  value  of  c*  using  N  =  3,  the 
agreement  with  observation,  though  good,  was  not  quite  as  satisfactory.  For 
this  reason  the  theoretical  prediction  was  probed  for  values  of  c*  =  0.9  and  1.0 
as  well. 

The  results  are  shown  in  the  following  Tables: 


c*  =  0.9 


uP,  U  in  mm. 


U  =  1 

U  =  2 

U  =  3 

0.010 

0.085 

0.219 

u| 

0.063 

0.293 

0.635 

0.280 

0.821 

1.463 

c"^  =  1.0 


U  =  1 

U  =  2 

U  =  3 

0.013 

0.104 

0.251 

uf 

0.073 

0.316 

0.677 

0.298 

0.863 

1.519 

3-7 


c*  =  1.2 


U  =  1 

U  =  2 

U  =  3  ' 

0.02 

0.142 

0.315 

uf 

0.094 

0.361 

0.760 

0.334 

0.947 

1.630 

Notation:  uf  =  uP  I  ^74 ,  uf  =  uP  I  ^72 ,  uf  =  uP  1 31,74. 


A  plot  of  uP  versus  x  for  various  values  of  U  is  shown  in  Figure  2, 
while  plots  of  uP  versus  U  at  different  points  along  the  height  of  the  cuboid 
are  shown  in  Figure  3.  Experimental  points  are  also  shown  at  x  =  0.56  h. 

An  Exact  Solution 


In  Figure  3,  the  decay  of  uP  with  respect  to  x  is  sufficiently  fast  to 
warrant  the  question  of  whether  the  fixed  end  of  the  cuboid  is  not  sufficiently 
far  from  the  free  end  to  warrant  the  treatment  of  the  length  of  the  cuboid  as 
"infinite".  We  show  in  this  sub-section  that  this  is  in  fact  the  case.  However, 
dz  is  still  set  equal  to  dU/h  where  h  is  a  characteristic  length,  the  actual  length 
of  the  cuboid.  The  value  of  h  is  relevant  only  insofar  as  it  concerns  the 
correlation  of  c*  to  the  value  of  the  diffusion  constant  a/h  in  the  generic 
equation  (3.1)  which  now  becomes: 


a  (9"uP  _  (9uP 

h  ^2 


(3.29) 


To  proceed  with  the  solution  of  this  equation  we  take  its  Laplace  transform  in 
light  of  the  initial  condition  uP(o)  =  0.  Thus: 


'^a^(9^u^  — p 
—  - =  p  u*^ 

IhJ 


(3.30) 


where  p  is  the  Laplace  transform  parameter. 


3-8 


If  c*  the  length  of  the  cuboid  is  infinite,  then  the  boundedness 
condition  on  u  as  x  — >  oo  reduces  the  solution  of  Eq.  (3.30)  to  the  form 


(3.31) 


Now  at  X  =  0,  u  =  U  and  thus  u  =  l/p^.  Hence, 


However,  we  note  that 


(3.32) 


(3.33) 


where  L'^  is  the  inverse  Laplace  Transform  operator  and  erfx  is  the  Error 
function.  Thus: 


Uf 


uP  = 


1 

ViF 


HU' 


J 


(3.34) 


which  is  a  closed  form  (exact)  solution  for  uP. 
Alternatively, 


duP 

"dlT 


’Vu 


(3.25) 
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Letting 


duP 

'W 


=  1  -  erf  k 


(3.26) 


Now  it  follows  from  Eq.  (3.28)  that 


\ 


=  1 


Hence, 


and  k  =  2. 


Thus, 


duP 

dU 


( 

1-erf 

V 


iy_ 

Vu 


\ 

y 


(3.27) 


(3.28) 


Below  is  a  table  of  values  for  duP/dU  for  both  the  finite  difference  solution 
when  N  =  3  and  the  exact  solution  given  by  Eq.  (3.27)  with  the  asymptotic 
assumption  of  infinite  length.  The  comparison  is  given  at  the  two  stations 
X  =  h/4  and  x  =  h/2.  In  the  case  of  x  =  h/2  a  comparison  with  the 
experimental  values  is  also  given.  The  agreement  is  close 

x  =  h/4 


U 

u'(N=3) 

u'(Exact) 


1 

0.47 

0.45 


2 

0.61 

0.61 


3 

0.69 

0.72 
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x=h/2 


u 

1 

2 

3 

u'(N=3) 

0.132 

0.32 

0.42 

u'(Exact) 

0.134 

0.32 

0.42 

u'  (Exp.) 

0.16 

0.33 

0.42 

duP 


Case  (2).  Unloading- 

One  set  of  experiments  by  Gill  [4]  probed  the  axial  displacement  field 
under  conditions  of  unloading.  The  result  of  one  typical  such  experiment  is 
shown  in  Figure  4.  The  essential  finding  is  that,  upon  unloading  the  relative 
uP  versus  U  is  essentially  a  straight  line  parallel  to  the  line  oB.  We  observe 
that  the  line  oB  is  the  response  appropriate  to  a  local  theory  ("local  behavior") 
which  in  the  context  of  the  present  theory  can  only  be  an  elastic  response  in 
the  light  of  Eq.  (3.39). 

The  inference  is  drawn  as  follows.  During  elastic  deformation 
duP/dz  =  0,  which  implies  (in  view  of  Eq.  (2.1))  that. 


n  =  0  (3.28) 

Hence,  during  unloading,  all  particles  are  immobile  and  the  deformation  of 
the  aggregate  results  from  the  elastic  expansion  of  the  grains  themselves. 


Thus  during  unloading  there  is  no  mobile  phase.  The  entire  aggregate 
domain  is  an  elastic  phase  and  hence  Qi  =  0  in  Eq.  (2.34).  One  concludes  that 
during  unloading  (where  uf  and  \\i  are  measured  relative  to  the  reference 
state  at  the  onset  of  unloading) 


=  0 


in  D 


(3.29) 
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In  the  one-dimensional  case, 


(7"  = 


^  dy/^ 
v^x  y 


=  constant  =  a 


(3.30) 


and  hence 


(9u 

=  —  =  cons  tan  t 
ox  ox 


(3.31) 


where  u  is  measured  relative  to  the  reference  state  at  the  onset  of  unloading. 
Thus  the  theory  is  local  and  u^  is  a  linear  function  of  x,  a  fact  that  agrees  with 
the  experimental  data  of  Gill.  See  Figure  5. 


Solution  of  Eq.  (3.31)  yields  the  result; 


u  =  uf 


V 


-1 

hy 


(3.31a) 


since  u  =  U  at  x  =  h. 

The  Stress  Response  Under  Uniaxial  Compression 
Loading  Case 

We  begin  with  Eqs.  (2.45)  and  (2.22b)  to  find  the  following  relation  for 
o; 

<7=-Aof(cT)s(w)^  (3.32) 

n  ox 

The  constants  (N/n)  and  Aq  in  conjunction  with  the  function  f  define  a  new 
function  F(a)  =  (N/n)Aof(a),  so  that  a  may  be  written  in  the  form 
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CT  =  F(<T)g(w)-^ 
d\ 


(3.32a) 


We  point  out  that  because  the  theory  is  non-local  w  varies  pointwise  along 
the  height  of  the  cuboid.  However,  because  of  the  small  strain  involved,  the 
porosity  does  not  vary  substantially  during  deformation.  Therefore  the 
dependence  of  g  on  w  is  weak.  For  this  reason  and  to  simplify  the  analysis  in 
this  investigative  stage  of  the  work,  we  regard  w  as  the  average  porosity  of  the 
cuboid. 

To  express  the  effect  of  w  analytically  we  reason  that  the  porosity 
depends  on  the  hydrostatic  strain.  This  is  true  in  principle  even  though  we 
realize  that  porosity,  strictly,  is  due  to  positional  arrangement  of  the  grains. 
Thus  we  set 

w  =  w(€)  (3.32b) 

and 

g[N(€)]  =  G(€)  (3.32c) 

Thus,  Eq.  (3.32)  becomes, 

a  =  F(CT)G(s)^ 

dx 

This  is  an  implicit  relation  between  a  and  3uP/3x. 

Remark.  There  exists  a  philosophical  question  of  determination  insofar  as 
the  function  G(e )  is  concerned.'  The  difficulty  lies  in  the  fact  that  the  strain  in 
the  reference  state  is  not  known  and  is  conveniently  assumed  to  be  zero. 
However,  following  a  history  of  loading  and  unloading  to  zero  stress,  there 
remains  a  residual  strain  which  requires  a  change  in  porosity.  Thus  during  a 
reloading  history  from  that  point,  €  cannot  be  set  equal  to  zero,  but  to  e  o 
which  denotes  the  residual  strain.  If  such  a  history  were  not  known  then  e  © 
would  be  an  unknown  constant.  Hence  G(€ )  is  in  fact , 
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G  —  G(e Q+  g)  . 


(3.33) 


where  €  o  is  a  function  of  previous  history  and  is  not  known.  The  problem  is 
obviated  if  G  is  such  a  function  that: 

G(€o  +  e)  =  kG(€)  (3.34) 

where  k  =  k(€  o)  and,  therefore,  a  constant.  Equation  (3.34)  is  satisfied  if 

G(€)  =  el5e  (3.35) 

where  P  is  a  material  constant. 

Thus  Eq.  (3.32)  becomes 

<7  =  F(  cr)e^^  (3.36) 

Since  the  constant  k  in  Eq.  (3.34)  can  be  absorbed  in  the  function  F,  one  can  set 
without  difficulty  €  o  =  0  in  the  reference  state,  i.e.,  prior  to  the 
commencement  of  a  strain  history.  We  shall  see  that  experiments  by  Gill  [6] 
do  in  fact  corroborate  the  form  of  G  depicted  in  Eq.  (3.35). 

The  Stress  at  the  Moving  Boundary 

The  experimental  results  by  Gill  [6]  were  presented  by  plotting  the  stress 
c  versus  the  displacement  U  at  the  moving  boundary,  but  also  versus  the 
"strain"  U/h.  Of  course  the  two  presentations  are  equivalent. 

A  plot  of  the  experimentally  obtained  axial  stress  plotted  versus  axial 
strain  is  shown  in  Figure  5,  for  a  history  of  three  loading  paths  oai,  a2a3,  a4a5 
and  three  unloading  paths  aia2,  a3a4,  a5a6-  In  our  discussion  of  the 
experimentally  obtained  mean  particle  displacement  in  the  neighborhood 
close  to  the  midpoint  of  the  cuboid  we  reasoned  that  unloading  behavior  is 
elastic. 
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This  conclusion  implies  that  the  unloading  paths  are  representative  of 
elastic  behavior.  If  that  is  the  case  then  the  curves  a2ai,  a4a3  and  a6a5  should 
all  be  part  of  one  and  the  same  curve.  This  is,  in  fact,  the  case.  The  elastic 
response  is  represented  quite  accurately  by  the  analytical  form 

/  p\6.25 

)  (3.37) 

as  shown  in  Figure  6.  The  analytical  representation  is  not  of  central 
importance.  Of  importance  is  the  fact  that  all  three  unloading  stress-strain 
curves  are  part  of  one  and  the  same  curve. 

We  observe  that,  in  a  manner  consistent  with  our  proposed  physics  of 
granular  motion,  loading  (reloading)  is  not  an  elastic  process  because  the 
mobility  is  not  zero.  Thus  the  curve  aoai  does  not  coincide  with  the 
unloading  curve  aia2  and  neither  do  the  other  corresponding  curves.  Clearly 
one  cannot  depict  such  results  by  means  of  a  yield  surface! 

To  determine  the  stress  during  loading  at  the  moving  boundary  of  the 
cuboid,  we  make  use  of  the  "infinite  length"  idealization,  which  as  we  have 
seen,  gave  an  accurate  displacement  distribution  in  that  it  matched  the  finite 
difference  solution  obtained  using  he  actual  length  of  the  cuboid. 

To  this  end  we  begin  with  Eq.  (3.32)  and  differentiate  both  sides  with 
respect  to  x  to  find 


(  \ 


Specifically  at  x  =  0,  i.e.,  the  moving  boundary. 


(3.38) 
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(3.39) 


Upon  taking  inverse  transform  of  Eq.  (3.39)  we  find: 


or 


duP  1 
dx 


■V^u 


duP 


dx 


•(€)^2 


(3.40) 


(3.41) 


where  e  =  U  /h. 

At  this  point  we  utilize  Eq.  (3.32)  and  set 


(3.42) 


to  obtain  the  following  expression  for  the  stress  at  the  free  boundary: 


<7  =  F*(<T)e^/*^ 


(3.44) 


Equation  (3.44)  is  the  basic  equation  that  will  be  used  for  the  theoretical 
depiction  of  the  loading  curve  aoai,  i.e.,  the  initial  loading  curve. 

Discussion.  Every  point  on  the  stress-strain  path  in  Figure  6  must  lie 
on  a  solution  curve  of  the  partial  differential  equation  (3.29).  The  unloading 
paths  satisfy  Eq.  (3.29)  identically  since  on  those  paths  uP  =  0.  The  loading 
paths  must  satisfy  Eq.  (3.29)  but  with  different  initial  conditions. 
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We  observe,  however,  that  Eq.  (3.29)  remains  invariant  with  the 
"translation  transformation". 


U'  =  U-Uo 


(3.45a) 


(3.45b) 


Thus  Eq.  (3.41),  i.e.,  the  solution  for  path  oai  applies  to  path  a2a3,  if  Uo  and  uP 
are  the  values  of  these  variables  at  the  point  a2.  Hence, 


(3.46) 


However,  insofar  as  the  function  G  in  Eq.  (3.32)  is  concerned,  6  and, 
therefore,  U  are  related  to  the  actual  total  change  in  the  mean  porosity  and 
therefore  must  be  measured  relative  to  the  initial  reference  point  0.  Thus  for 
the  loading  paths  a2a3  and  a4a5,  Eq.  (3.44)  will  read 

<T  =  F*(c7)e^^(e')^  (3.47) 


there  €  is  the  cumulative  strain  (i.e.,  the  strain  measured  from  the  reference 
point  0)  and  e '  is  measured  from  the  point  of  reloading. 

The  comparison  between  theory  and  experiment  is  shown  in  Figure  7 
with  P  =  0.06  and  F*(a)  as  shown  in  Figure  8.  How  P  and  F*  are  determined 
will  be  the  subject  of  the  next  subsection. 
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Plastic  Displacement  (mm) 


Fig.  2.  Theoretical  Plastic  Displacement  vs  Height 
for  Various  Values  of  Surface  Displacement 


Plastic  Displacement  (mm) 


1.5 


•  Experiment 
—  Theory 


3: 


Fig. 3.  Plastic  Displacement  vs  Surface  Displacement 
at  Various  Height  Locations 
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Displacement  of  Loading  Cap  (m) 


Macroscopic  Strain 

Displacment  of  the  Average  Particle  Centroid; 

A;  ADJACENT  TO  THE  LOADING  CAP. 

B;  AT  0.56h  if  a  linear  displacement  distribution 
IS  ASSUMED  (local  THEORY) . 


Fig.  4.  Experimental  Displacement  Points  Near  the  Mid-Point 
and  Comparison  With  the  Theoretical  Line 
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Fig.  5.  Experimentally  Obtained 
Stress-Strain  Response 
Under  Axial  Strain  Conditions 
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Fig.  7.  Experimental  and  Theoretical 
Stress-Strain  Response 
Under  Axial  Strain  Conditions 
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THE  FUNCTIONS  F  (a)  and  H(a) 


10  20  30  40  50 

STRESS  (MPA) 


Fig.  8.  The  Functions  F*  and  H  vs  Stress 
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SECTION  4 


DETERMINATION  OF  MATERIAL  CONSTANTS 


ONE-DIMENSIONAL  HISTORIES 

We  begin  this  Section  with  the  method  of  determination  of  material 
constants,  and  functions  in  one-dimensional  strain  histories. 


Case  (il  Loading  Histories 


The  relevant  equation  is  Eq.  (3.29)  —  which  determines  the 
displacement  field  in  the  cuboid  —  and  its  companion  equation  (3.1a),  which 
we  repeat  here  as  Eqs.  (4.1  and  (4.1a): 


a  (9^uP  _ 
h  (9x2  "lu 


(4.1) 


a  = 


^  A  ^ 
V  J 


(4.1a) 


and  the  equation  for  the  determination  of  the  stress,  i.e., 
<T  =  F*((7)e^^(€')^ 


(4.2) 


which  is  Eq.  (3.47).  We  note  that  in  the  case  of  loading  histories 
uP  =  u  ,  eP=€ 


(4.3a,b) 


Case  (ii).  Unloading  Histories 

The  relevant  equations  here  are: 
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(4.4) 


<T=cr(e^) 

Method  of  Determination  of  the  Constants  and  Functions 


(4.5) 


Case  (i).  Loading  Histories 

The  constants  and  functions  to  be  determined  are  the  constant  "a"  in 
Eq.  (4.1),  the  constant  P  in  Eq.  (4.2)  an  the  function  F*(a)  in  Eq.  (4.2). 

The  Constant  "a".  Given  a  granular  medium  in  the  form  of  a  cuboid  it 
w'ill  not  be  known,  ab  initio^  that  the  solution  of  Eq.  (4.1)  can  be  obtained  by  an 
infinite  length  approximation,  because  this  will  depend  on  the  actual  height 
h  of  the  cuboid,  and  the  friction  coefficient  bo  as  well  as  the  elastic  constant 
Ao,  in  accordance  with  Eq  (4.1a).  Thus  the  closed  form  solution  given  in 
Eq.  (3.25)  cannot  be  used. 

The  approach,  therefore,  is  to  solve  Eq.  (4.1)  by  a  finite  difference 
technique  with  a  three  node  approximation,  one  node  at  x  =  0,  one  at  x  =  h/2 
and  one  at  x  =  h.  The  boundary  conditions  are  such  that  uP(o)  -  0,  uP(h/2)  = 
uP  and  uP(h)  =  U.  The  solution  is  given  by  Eq.  (3.27),  i.e.. 


where  c  =  8a /h^. 

Thus  if  uP  is  measured  at  or  close  to  the  midpoint  of  the  cuboid  the 
constant  c  can  be  determined  by  matching  the  right-hand  side  of  Eq.  (4.6)  to 
the  measured  relation  between  uP  and  U  as  in  Figure  9.  The  value  of  c  thus 
found  as  0.6. 
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With  this  value  of  c  at  hand,  Eq.  (4.1)  was  then  solved  more  accurately 
using  five  nodes  in  the  finite  difference  scheme,  as  discussed  in  Section  3. 
The  solution  in  that  section  was  given  in  terms  of  values  of  c*(=  2  c)  equal  to 
0.9,  1  and  1.2,  and  was  determined  that  c*  =  1  was  the  optimal  value  of  c*  (i.e., 
c  =  0.5)  for  the  most  accurate  theoretical  representation  of  the  observed  data. 
Since  c*  =  16  a/h^  the  constant  "a"  was  thus  determined. 

The  Constant  p 

In  Figure  10  we  highlight  the  loading  stress-strain  curves  oai,  a2a3,  a4a5. 
The  horizontal  line  ABC  is  a  constant  stress  line  that  intersects  the  curves  oa, 
at  A,  a2a3  at  B  and  a4a5  at  C.  Thus  ga  =  ob  =  oc. 

In  reference  to  Eq.  (4.2)  we  note  the  following  value  of  e  c  in  units  of 
1  percent  (0.01): 

€(A)  =  2.2  ,  ,e(B)  =  3.0  ,  g(C)  =  (5.4). 

We  also  note  that 

e'(A)  =  2.2  ,  ,e'(B)  =  3.0  ,  €'(C)  =  (1.6). 


In  view  of  Eq.  (4.2),  we  have  the  following  relations: 


(4.7a) 

GB  =  F*{GB)e^^<^^(€fe)>^ 

(4.7b) 

(4.7c) 

Thus  it  follows  that,  since  ga  =  ob  =  oc, 
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(4.8a, b,c) 


Hence  typically 


p{e  (B)-  G  (A))  = 


(4.9) 


Six  values  of  P  for  six  such  computations  in  accordance  with  the  above 
procedure  are  shown: 


P  =  {0.046, 0.050, 0.060, 0.060, 0063, 0.066} 


An  average  P  was  then  set  at  0.060. 
The  Function  F*(a) 


The  function  was  determined  directly  for  an  Eq.  (4.2)  and  the  loading 
curve  a4a5  of  Figure  10.  Evidently, 

F"(a)  =  f7e‘'^^(€')^  (4.10) 

where  €  and  e  ’  are  measured  in  units  of  1  percent. 

Since  every  term  on  the  right-hand  side  of  Eq.  (4.10)  is  known  at  every 
point  on  the  curve  a4a5,  then  F*(a)  can  be  determined  for  all  o  on  that  curve. 
The  function  F’^(a)  so  determined  is  shown  in  Figure  10.  Also  shown  in  that 
Figure  is  the  function  H(a)  where 

H(<t)  =  ^:^  (4.11) 

F^(ct) 


Hence, 


and  therefore. 


(4.13) 


The  importance  of  Eq.  (4.13)  is  that  it  is  an  explicit  equation  for  the  stress  for 
all  loading  curves,  emanating  at  a  =  0.  In  fact  this  equation  was  used  to 
predict  the  stress  response  for  the  curx^es  oa,  and  a2a3,  a  result  already  shown 
in  Figure  8. 

Determination  of  Material  Constants  ad  Function  in  Three  Dimensions 
1.  The  Mobile  Phase 


The  three-dimensional  motion  of  the  mobile  phase  is  given  by  Eqs. 
(2.54)  and  (2.57),  i.e.. 


11  u-,i|  =  b 


dz 


(4.14) 


and 


(A+2/y)(/),ii  =  b^  (4.15) 

Recall  that  the  displacement  field  Ui  is  given  by  Eq.  (4.16): 

uP  =u[-i-ui  (4.16) 

where  uf  is  an  isochronic  field,  i.e.,  uf^  =  0  while  Ui  is  irrotational,  i.e.,  Ui,j  = 
Uj  i,  and: 

u,  =  (t>,i  (4.17) 
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The  Irrotational  Field 


In  the  uniaxial  experiment  discussed  in  Section  3, 
u^  =  uf’(xi) ,  u^  =  0 ,  u^  =  0  . 

Thus 

Uj  ^  =  U2  2 ,  Uj  ^  =  U2  3  =  Uj  2  • 

and  the  field  is  irrotational. 

In  this,  uniaxial,  case  Eq.  (2.55)  may  be  applied  directly.  Thus 


(A  -h  2/y)Uxx  t)-— 

OT 

(4.18) 

where  we  have  set  xi  =  x.  Thus,  in  Eq.  (2.41) 

A  =  A  +  2jL/ 

(4.19) 

and 

X+lfd  - 

a-  b 

(4.20) 

as  in  Eqs.  (3.32d)  and  (3.35). 

Also, 

A  =  A„F(<7)e^® 

(4.21) 

b  =  b<,F((j)e^® 

(4.22) 
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Thus  the  uniaxial  test  determines  the  constants  and  functions  of  the 
irrotational  motion,  as  discussed  earlier  in  this  section. 

The  Isochronic  Field 


Here  the  appropriate  experiment  is  a  shear  experiment  such  as  the  one 
conducted  in  a  "sheer  box".  See  Ref.,  7  for  instance.  Of  note  is  the  fact  that  a  is 
common  to  both  the  isochronic  and  irrotational  fields.  Thus  the  function 
F(a)  and  the  constant  p  are  determined  from  the  irrotational  test.  There 
remains  the  constant  a*  where 


(4.23) 


and  governs  the  isochronic  motion  given  by  Eq.  (4.14). 

The  procedure  for  determining  a*  is  precisely  the  same  as  that  for 
determining  a(a)  in  Section  3,  by  measuring  the  displacement  at  the  half¬ 
height  of  the  shear  box  as  a  function  of  the  surface  displacement  and  then 
using  Eq.  (4.6). 

THE  ELASTIC  PHASE 


We  proceed  on  the  assumption,  normally  made  in  soil  mechanics,  that 
the  elastic  shear  response  is  linear  while  the  hydrostatic  response  is  non¬ 
linear.  Thus 

Sij=2/ye^  (4.24) 

cr=K(G^)  (4.25) 

where  K  is  a  function  of  the  elastic  hydrostatic  strain. 
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The  Shear  Modulus  |i 


Under  pure  sheer  unloading  conditions  let  s  be  the  shear  stress  and  e^ 
the  corresponding  elastic  shear  strain.  Then 

2iJ=s/e^  (4.26) 

and  hence  p  is  found. 

The  Function  K(g 

This  may  be  found  by  means  of  the  uniaxial  strain  test  under 
unloading  conditions.  We  note  that,  in  view  of  Eq.  (4.24) 

aj  -  (72  =  2/y  €]  (4.27) 


since  €2  =  0.  Also  in  view  of  Eq.  (4.25) 


^€i^ 


(ai+2cT2)  =  3K  ^  =K’^(€i) 


Eliminating  02  from  Eqs.  (4.27)  and  (4.28); 

=  +K=^(ei) 

Hence  K*  and  thus  K(e )  may  thus  be  found. 


(4.28) 


(4.29) 
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Plastic  Displacement  (mm) 


Fig.  9.  Plastic  Displacement  vs  Surface  Displacement 
at  Mid-Point 
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SECTION  5 


DISCUSSION 

There  is  ample  evidence,  both  experimental  [6]  and  theoretical,  (as 
discussed  by  the  author  in  past  work  [1]  and  in  the  current  study  and  by  other 
authors  [7])  to  the  effect  that  the  mechanical  response  of  granular  domains, 
under  surface  tractions,  is  non-local.  The  precise  nature  of  the  non-local 
behavior  is  a  matter  for  further  research. 

In  this  study  and  in  the  work  that  preceded  it  [1-3]  we  have  identified 
three  causes  for  non-locality:  (a)  internal  friction  giving  rise  to  non-local 
thermodynamic  formulations,  (b)  configurational  entropy  and  temperature 
giving  rise  to  diffusive  motion  of  the  "thermal"  kind  and  (c)  mobility  which 
is  most  likely  to  be  a  field  with  its  own  field  equation. 

RECOMMENDATIONS 

Mindful  of  the  need  to  lend  practical  value  to  constitutive  equations, 
further  research  is  needed  to  identify  the  practical  circumstances  where  a  local 
theory  will  suffice,  but  also  reveal  "pathological  circumstances",  such  as- 
banding  and  instability,  where  clearly  a  non-local  theory  is  needed  for  the 
physical  understanding  and  prediction  of  such  pathological  events. 

Continued  research  in  the  study  of  non-local  effects  in  the  motion  of 
granular  domains  is  therefore  strongly  recommended. 
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